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SADDLEPOINT  APPROXIMATIONS  FOR  THE  COMBINED  PROBABILITY  AND 
JOINT  PROBABILITY  DENSITY  FUNCTION  OF  SELECTED  ORDER  STATISTICS 

AND  THE  SUM  OF  THE  REMAINDER 


INTRODUCTION 


Detection  and  location  of  weak  signals  in  random  noise  are  frequently  accomplished  by  the 
ordering  of  the  random  variables  (RVs)  in  a  measured  dataset,  followed  by  an  investigation  of 
the  locations  and  statistics  of  several  of  the  largest  RVs  under  consideration.  Also  of  interest  are 
the  remaining  smaller  RVs  in  the  dataset,  which  can  be  used  to  estimate  the  background  noise 
level  and  to  form  a  basis  for  normalization,  thereby  realizing  a  constant  false  alarm  processor. 

In  this  report,  the  original  dataset  {x„ }  is  composed  of  N  independent  RVs  with  arbitrary 
probability  density  functions  (PDFs)  {p„(x)}.  This  dataset  is  ordered  into  another  dataset  of 

dependent  RVs,  each  with  a  different  PDF.  From  this  latter  set,  the  M-  1  largest  RVs  are 
selected.  The  sum  of  the  remaining  N+  1  -  M  RVs  is  then  computed,  giving  a  total  of  M 
dependent  RVs.  The  joint  M-th  order  PDF  of  these  M  dependent  RVs  is  one  of  the  quantities  of 
interest. 

For  convenience,  in  this  report,  the  largest  RV  in  set{x„}  is  denoted  by  z,,the  second- 
largest  byz2, ...  ,the  M- 1  largest  byzM_15  and  the  sum  of  the  remaining  RVs  by  zM .  Thus,  the 
first  M-  1  RVs  satisfy  the  restrictions  that  z,  >  z2  >  •  •  •  >  zM_, ,  while  the  last  RV  satisfies  the 
restriction  that  zM  <  (N  + 1  -  M)  zM_, .  In  two  earlier  studies,  detailed  in  references  1  and  2,  the 
combined  probability  and  joint  PDF  of  these  M  RVs  was  derived  in  a  form  involving  a  one¬ 
dimensional  Bromwich  contour  integral  in  the  moment-generating  function  (MGF)  domain. 
Here,  a  saddlepoint  approximation  (SPA)  to  that  contour  integral  is  derived  and  applied  to 
several  typical  sets  of  RVs.  In  addition,  the  first-order  correction  term  to  the  SPA  is  developed. 
This  development  is  based  heavily  on  the  results  and  notations  of  references  1  and  2. 


SADDLEPOINT  APPROXIMATION  FOR  COMBINED  PROBABILITY  AND  JOINT 

PROBABILITY  DENSITY  FUNCTION 


The  combined  probability  and  joint  PDF  of  interest  is  given  in  reference  2,  equation  (51). 
namely, 

M-2  M- 1 

(^1  *  •  •  •  •>  >  •  •  • » =  nw.-ui  n<A.,<*>» 

m=l  j= i 

1  f  /M-\ 

x  -j2a  rA  *xp(~A  Zm  }  A )/  n  ,*)} ,  (1) 

C  /  7=1 

where  (/( )  is  the  unit-step  function,  C  is  the  contour  in  the  MGF  domain  X,  P(z ,  A)  is  the 
product  function  defined  in  equation  (3)  of  reference  2,  and  the  c-function  is 

v 

c„  («,  A)  =  Jr/v  exp(Ax)  pn  (x)  for  n  =  1 :  N.  (2) 

-or 

Evaluation  of  the  top  line  of  equation  (1)  is  trivial;  however,  the  bottom  line  is  more  involved. 

An  alternative  version  of  the  bottom  line  of  equation  (1)  is 

~  \dX  exp(-/l  zM)  Y\{cn{zM„X)},  (3) 


where  the  notation  neX denotes  n  =  1 : N  except  n  *  m1 ,  m2 , . . . ,  mM . 
The  logarithm  of  the  integrand  in  equation  (3)  is  defined  as 


A(X)  =  log  exp(-A  zm)Y[{c„{zm  „/l)} 

ncX 

=  ~A  +'Zl°glc„(=M  nA)l 


The  first  two  derivatives  with  respect  to ,4  follow  as 

A\X)  =  -z.,  +V  C-^Z"  '-1^1 
tyCn{- A,  rX) 


*^2  )  =  y  ^  ^  ( -A/  pj_)  [C,|(~A? 
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The  saddlepoint  (SP)  Xs  of  A.(A)  (and  the  integrand  of  equation  (3))  is  located  where 
A  (As )  =  0,  As  =  As  (zM ,  ), 

which  is  independent  of  variables  {zm},  m  <  M-  1,  but  which  does  depend  on  integers  {mj}, 
j  -  \  :M-  1.  The  standard  SPA  to  equation  (3)  is  given  by 


expM,  zM 


SPAO 


neX _ 


[2^AU>] 


exp[A(/ls)] 

P^AU)]1'2 


The  first-order  correction  to  the  SPA  is 


SPALSPAO  >  + 

8  A'U)2  24  A'(4)! 


and  requires  two  additional  derivatives  beyond  those  given  in  equations  (5)  and  (6). 


EXAMPLE  1:  EXPONENTIAL  RANDOM  VARIABLES 

From  pages  A-l  and  A-3  of  reference  1,  there  follows  the  individual  PDF 


1  f  X  | 

p„(x)  =  — exp - forx>0,  n  =  l:N 

K  b„ 


and  the  corresponding  o function 

l-exp(-u/b„ +uA)  ,  u 

c„(u,A)  = - - for «  >  0;  c„(w,l/6„)  =  — . 

1  -b„  A  bn 


Random  variable  x„  has  mean  b„  for  n  =  1  :N. 


At  this  point,  it  is  useful  to  define  function 


L(z)  =  log  f 1  CXp(  for  all  -  *  0;  L( 0)  =  0. 


Then,  using  equation  (11), 


\og[cn(u,A)]  =  \og  f  Ul  ^-uA  , 
{ b.,  b„ 


and  equation  (4)  yields 


(14) 


A(A)  -  -XzM 


Z  log 


N\ 

lm- i 


mX 


v  K  j 


+  L 


f  . 

*A/-1  _  „ 

U  LM-' 
V  *» 


There  follows. 


A'W  =  -z„-r„.,  £i' 


ncX 


cM-\ 

v  K 


‘■M-l 


X 


A,*V)=(--v-,)‘Ei< 


(*) 


*A/-1 


^A/-l  ^ 


for  A:  >  2. 


(15) 


The  derivatives  required  of  equation  (12)  use  equations  (7),  (8),  (9),  and  (15),  and  are  given 
by 


L'{z)  =  -I  +  — L-  £  =  exp(z), 
z  E  —  1 


2  £(£  +  1)  (,6) 

*3  +  (£-l)3’ 

6  £(£2+4£  +  l) 

~  z4  (£-l)4 


The  values  of  the  functions  above  at  r  =  0  are  to  be  taken  as  the  limits  as  s  approaches  zero.  All 
four  of  these  functions  are  analytic  functions  of  z  on  the  real  axis.  However,  great  care  must  be 
taken  in  dealing  numerically  with  the  removable  singularities  at  z  =  0,  especially  for  the  higher 
derivatives.  These  issues  are  discussed  in  appendix  A.  The  high-order  poles  at  -  =  i2nk, 

k  =  ±1,  ±  2,. . do  not  interfere  with  the  search  for  the  SP  Xs  of  equation  (7)  on  the  real  A-axis. 

That  is,  A'(X)  in  equation  (15)  encounters  purely  real  arguments  for£'(^)in  equation  (16)  when 
X  is  real. 
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EXAMPLE  2:  SQUARED-GAUSSIAN  RANDOM  VARIABLES 
From  pages  A-2  and  A-3  of  reference  1,  the  individual  PDF  is 


Pr,(x)  =  — - TTT7T  for  x  >  0 

(2  7rxbn) 


and  the  corresponding  c-function  is 


c„(u,A)  = 


1-2  U 


for  u  >  0;  c„  u, 


Random  variable  x„  has  mean  b„  for  n  =  1  :N. 


At  this  point,  the  relevant  function  is  defined  as 

L(z)  =  log  Cr  -  forallz^O;  Z,(0)  =  log(2/V^r), 

Vz 


which  is  entire.  Then,  using  equation  (18), 


log[c„  (u,  A)]  =  i  logf  -~T~\  +  L\~ —  u A 1 
2  l  2b„  l  2 j 


and  equation  (4)  yields 


AW  =  -^-«+^IlogfeL  +1^  ~ 

z  ntX  \LOn  )  nsX  V  Z 


There  follows 


A  (/l)  -  z^  -  zM_,  ^  I  ZM-1  ^  ’ 


'  2* 


A“'w- Z1'1'  for*22- 


The  derivatives  required  of  equation  (19)  use  equations  (7),  (8),  (9),  and  (22),  and  are  given 


l'V)=~  +  G,  C»  «*-»>  , 

2*  ybrz  erf (slz) 


1  i  o  5^ 

ifw= A- 1+—  g-g2, 

zz  v  -  y 

^(-)="V  +  fl  +  -  +  ^V  +  f3  +  — )g2+2G3, 


,  3  f  1.5  2.25  1.875V  V  7  3.75^  ,  f  6^  , 

1  (^)  =  7T“  1  +  —  +  ^-  +  — —  G-  7  +  — +  — —  G1  -  12  +  -  G3  -6G  . 
z  \  z  z  z  )  \  z  z  )  \  z ) 

The  values  of  the  functions  above  at  z  =  0  are  to  be  taken  as  the  limits  as  z  approaches  zero.  All 
four  of  these  functions  are  entire  functions  of  z.  However,  great  care  must  be  taken  in  dealing 
numerically  with  the  removable  singularities  at  z  =  0,  especially  for  the  higher  derivatives. 
These  issues  are  taken  up  in  appendix  B. 


EXAMPLE  3:  GAUSSIAN  RANDOM  VARIABLES 

From  pages  A-l  and  A-3  of  reference  1,  the  individual  PDF  is 


,  ,  1  If  x -a 

pAx)=^'*>  Vv 


and  the  corresponding  c-function  is 


h  H  b 


c„  (u,  A)  =  exp]  an  A +  \  b2  A2 1 of  — ^  -  b„  A 1 

V  2  J  {  K 

where  O  is  the  cumulative  distribution  function  of  a  normalized  Gaussian  RV.  That  is, 


0(0  =  jdu0(u)=  jdu  (2nyV7  exp(-u2/l). 


Random  variable  x„  has  mean  an  and  standard  deviation  bn  for  n  —  1  ‘.N  (see  equation  (24)). 
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At  this  point,  the  relevant  functions  are  defined  as 


1  1.2  l2  ,  At  5\\  At  1\  _  ^  ^ 


1(A)  =  aA  +  -bzA2  +  log  ®(A(A)),  4f(A)  = - b  A,  A'(A)  =  -b, 

2  b 


which  are  entire.  Then,  the  first  derivative  is  given  by 


L'(A)  =  a  +  b2A-br(A), 


where 


rW"0uajj’  '■'W-M-O^W+rW]-  P») 

The  remaining  higher  derivatives  follow  upon  reapplication  of  equation  (29)  to  equation  (28), 
namely, 

L\A )  =  b2  -  b2  r(A)  [A(A)  +  r(A)], 

Lm(A)  =  b 3  r(A)  [1  -  J 2  (A)  -  3 ^(A) r(A)  -  2  r 2  (A)],  (30) 

L”\A)  =  6 4  r(A)  [3 ^(A)  -  A3 (A)  +  (4  -  7 A2 (A)) r(A )  - 12  A(A)r2(A)  -  6 r 3 (A)]. 

The  quantities  .4(A)  and  r(A)  are  given  in  equations  (27)  and  (29),  respectively. 

The  corresponding  derivatives  of  equation  (4)  are  now  available  as 

AW  =  -Xzu  +2  [a,  1  +  i*2  l2  +  logOW)],  4  =  A,  W  A, 


AW  =  -zM  +  z  K  +  bl  A-  -  K  rn  l  r„  =  r„  (A)  = 

nsX 

A'W  =  Z*.!-ZA»'iK  +  '-J. 

A'U)  =  24  'i  p  -  4  -34.  -2 


4>M„W) 


A "(2)  =  24  ^  PA  -  A3  +  (4  -  743)  r,  -12 A,  r;  -  6r.J], 
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EXAMPLE  4:  CHI-SQUARED  RANDOM  VARIABLES  OF  2k+2  DEGREES  OF 
FREEDOM 

From  pages  A-2  and  A-3  of  reference  1,  the  individual  PDF  is 

,  .  xk  exp(-x/b)  _  „  ,  . 

P„(x)  = - TjTi  -  for  x  >  0,  k  integer  >  0 

UfJ 

and  the  corresponding  c-function  is 

k  1 

l-exp  (ru/b„  +  uX)YJ-(ulb„-uZ)i 

cMX)=  (HT -  for">0; 


Random  variable  x„  has  mean  (k+\)  b„  for  n  =  1  :N. 
Define  function 


l-exp  (s)£zJ/f. 


L(z)  =  log 


for  r  *  0;  L( 0)  =  -  log((£  + 1)|) 


Then,  using  equation  (33), 

log  [c„(w,  Z)]  =  (k  + 1)  logf  ^  +l[--uZ 

\°J  \K 

and  equation  (4)  yields 


\(Z)  =  -ZzM  +  (*  +  !)£  log 


n  J  ncX 


There  follows 
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\^L'  *f  1  zM_xX  , 
heT  V  J 


A(^(2)  =  (-za/_1)^  forp>2. 

«£  v  K  J 


At  this  point,  it  is  necessary  to  consider  specific  cases  for  integer  k  in  equations  (32)  through 
(34).  The  case  of  h  =  0  was  discussed  in  equation  (10)  and  the  text  that  followed.  The  first  case 
of  interest  here  is  k=  1,  namely,  2k+2  =  4  degrees  of  freedom  (DOF).  The  function  of  interest 
and  its  first  four  derivatives  are 

L(z)  =  logf  -  ~  eXP(~f-)-^-+-!)"|  for  z*0,  L( 0)  =  -log(2), 

V  z  J 

L'(z)---  +  J  ,  £  =  ex p(z), 

4  E1  (-2  +  z)  +  E(4-z  +  z2)-2 
{Z>~  z 3  (E-l-zf 

s_  12  £3  (3-z)  +  4£2(-3  +  2z-z2)  +  £  (15-7z  +  z2  -z3)-6 

W'r4+  (E-l-zY 

The  second  case  is  k  =  2,  which  corresponds  to  six  DOF: 


l-exp(-z)  (l  +  z  +  z2/2)  )  .  fC. 

L{z )  =  log  - - - -  for  z  *  0,  £(0)  =  -log(6), 

V  2  J 

L\z)  =  E  s  exp(z),  Dsf-(l  +  z  +  z2/2), 

r(-)  =  J_+-^zi)zl=£ 

z2  '  2D2 

6  |  £2(4-8z  +  2z2)  +  £(-8  +  8z  +  8z2-2z3+z4)  +  4-6z2 -2z3 

r"(z)  =  ^-  +  — !— (4£3(-6  +  6z-z2)  +  8£2(9-8z2  +4z3  -z4) 
z4  8D4  V 

+  £(-72-  72r  +  68r2  +  16z3  -  14z4  +  2z5  - z6)  +  6(4  + 8z  - 4z3  -  z4 )) 
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The  third  case  is  k  —  3,  which  corresponds  to  8  DOF: 


L(z)  =  log 


1  ~  exp(-r)  (l  +  z  +  z2/2  +  z3/6) 


for  z  *  0,  1(0)  =  -log(24), 


4  z3 


I'(r)= ,  E  =  exp(z),  D=  £-(l  +  z  +  z2/2  +  z3/6), 


i'W  ■ -  p~ +  (£  (6 -  2.-)  -  (6  +  42  +  z!  )} 

i"(‘)  =  -4  +  ^7  («£’  (6  -62  +  2s)  +  £(-72  +  24r2  +  122s  -  324  +  2s) 

2  3o  U 

+  36  +  36z  +  6z2  -12z3  -6z4  - z5), 

(36£3  (6  - 1 8z  +  9z2  -  z3)  +  24 E2  (-27  +  54z  +  27 z2  -  15z4  +  6z5  -  z6) 

+  £(648-648z-1620z2  -1 188z3  +126z4  +162z5  +3z6  -21z7  +3z8  -z9) 

-  3(72  - 216z2  -  336z3  -  186z4  -  24z5  +  20z6  +  8z7  +  z* ))+ 

2 

The  method  of  handling  the  high-order  removable  singularities  of  equations  (38),  (39),  and 
(40)  at  z  =  0  is  discussed  in  appendix  C. 


SEARCHING  FOR  THE  SADDLEPOINT 

The  SP  location  of  the  integrand  of  equation  (3)  is  on  the  real  X  axis.  Furthermore,  the 
integrand  has  a  single  minimum  on  the  real  axis,  in  its  region  of  analyticity,  which  is  the  SP 
location  Xs.  This  will  be  demonstrated  for  arbitrary  PDFs  {pn(x)}  and  arbitrary  M.  It  is 
presumed  that  zm<  (N+\-M)  zm-  \  so  that  the  joint  PDF  of  interest  is  nonzero. 

Equation  (6)  is  the  starting  point  for  the  second  derivative  of  the  logarithm  of  the  integrand 
of  equation  (3).  In  general, 


u 

c{u,A)  =  Jdv  p(x)  exp(i  x), 

-or 

1/ 

c'(u.A)  =  p(x)  exp(,l  x)  x, 

-  or 

U 

c\u. A)  =  Jrfv  p(x)  exp(Ax)  x2 . 
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Then,  the  numerator  of  equation  (6)  can  be  expressed  as 


l2 


U  U  ll 

^dx  p(x)  exp(A  x)  jdx  p(x)  exp(Ax)  x2  -  ^dx  p(x)  exp(Ax)  x 


(42) 


But,  Schwartz’s  inequality  states  that 


[i 


dx  f{x)  g(x) 


-\2  r 
< 


fdxf2(x)  J  dx  g2(x) 


(43) 


with  equality  if  and  only  if/(*)  and  g(x)  are  proportional  for  all  x.  Upon  identification  of 


/ (x)  =  -J p(x)  exp(T  xf  2)  and  g(x)  =  yfpix)  exp(/l  x/2)  x. 


(44) 


it  follows  that  equation  (42)  is  positive  for  all  X.  Therefore,  every  numerator  of  equation  (6)  is 
positive,  meaning  that  A  "(A)  >  0  for  all  real  X.  That  is,  A(/L)  is  bowl-like,  for  real  A,  with  a  single 
minimum  in  its  original  region  of  definition  (convergence  of  its  integral).  Therefore,  the 
integrand  of  equation  (3),  namely  exp[A(A)] ,  has  the  same  property. 

The  above  conclusion  holds  for  all  values  of  zm  and  zm-i-  However,  the  actual  location  of  the 
minimum  of  A(A)  depends  on  both  of  these  variables,  as  can  be  seen  by  setting  equation  (5)  to 
zero  (see  equation  (7)).  Variables  {zm}  for  m  <  M-  1  have  no  effect  on  SP  location  A*. 
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RELATIVE  ACCURACY  OF  THE  ESTIMATE 
OF  THE  SADDLEPOINT  APPROXIMATION 


The  MGF  corresponding  to  PDF  p(u )  is  given  by 
p{A)  =  jdu  exp  (Au)  p(u ) 


and  the  inversion  formula  is 


p(u)  =  fdA  exp (-Au)  p(A), 

1  /.IT  * 


where  Bromwich  contour  C  is  a  vertical  line  from  -ioo  to  +ioo  in  the  region  of  definition 
(analyticity)  of  MGF  p(A),  as  obtained  from  convergence  of  the  integral  in  equation  (45).  The 
logarithm  of  the  integrand  of  equation  (46)  is  defined  as 

A(A,u)  =  -Au  +  x(A),  where  j(A)  =  log  p(A).  (4 

The  SPA  to  exact  PDF  p(u)  in  equation  (46)  is  then 


Ps(*)  = 


exp(-/^  u )  pjAJ 

V*X'(AS)\V2  ' 


where  SP  location  As  =  As(u)  satisfies  the  SP  equation,  z'(As)  =  u. 

Now,  in  the  neighborhood  of  the  actual  SP, 

Z'(A)  =  x'(As )  +  z”(As )  (A  -  As )  for  A  near  As , 
z'(A)-u  =  x'(As)(A-As). 

Let  a  relative  error  tolerance  Ra  be  imposed  on  the  accuracy  of  the  SP  location,  namely, 


X'(A)-u 


<  Ra;  «*0. 


Then,  the  solution  A  to  the  SP  equation  will  satisfy 

AH. 

v  x  (K ) 
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Therefore,  let  solution  Ac  (the  axis  crossing  actually  used  for  the  Bromwich  contour  C)  be 

I  I  \U\ 

A=Z  +  e,  where  Is]  <  — . 

II  zW 

Adopt  the  estimate  of  the  SPA  as 
_ /.a  exp (~Acu)n(Ac) 


Pc(U)  = 


\2n  X\K)\in 


in  analogy  with  equation  (48).  Then, 

exp(-Ac  u )  =  exp(-(As  +s)u)  =  exp(-As  u)  exp(-£  u )  =  exp(-li  u)(\-£  u), 


ju(Ac)  =  ju(As  +e)  =  m(As)  +  £  M\XS)  =  p(As )  1  +  £ 


=  )  [1 + £  zXK )] = P(K  )  0 + £ u). 


pW 

pas) 


exp(-Ac  u )  p(Ac)  =  exp (-As  u)  ju(As) (l-£  u  ), 

which  is  quadratic  in  £.  This  is  consistent  with  As  being  a  SP  of  the  integrand  of  equation  (46). 
Also, 

zW = z\K + £) = z"W+ £  zmW, 


r  3  vil/2  ^  r  n  /  o  yil/2  i  I  £  Z  (^5 ) 

L x  (^c)]  —  \x  )j  *  +  ~  "s  i  \ 

^  X  \Xsj 


Combining  these  results  with  equations  (53)  through  (55),  there  follows 
PcO*)  =  P,(M)  ]~f  X  to  order  £. 

L  2  z  ajj 

Therefore,  using  equation  (52),  the  relative  error  in  the  estimated  SPA  pdu)  is 


*  raj  < 

„  *.  « 

zm(A) 

2  Z'W,)  ' 

2 

zW! 
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Recall  that  SP  location  As  —  as(u)  is  a  function  of  u.  Therefore,  R\  can  have  a  complicated 
dependence  on  field  point  u.  The  third  derivative  of  x  at  As  only  needs  to  be  evaluated  once ,  after 
the  SP  location  has  been  determined.  Ratio  R\!Ra  is  called  the  error  factor. 

The  modified  correction  term  (CT)  to  the  SPA  is 
cr_l  5  X"(AC)\  ,  m 

(58) 

Upon  expanding  these  functions  according  to 

)  =  xXK  +  e)  s  z\*.)  +  e x"(As )  =  Z2+SX3, 

(59) 

Xm(K)  =  X3+EX^  X”XK)  =  X<+£X^ 
there  follows 


CTz-£ f-  !  +  £■[  — -2 

8  x2  u4 


-AZl  i  +  £  2Il-3^-|. 

Xt  X2  J J  24  x2  {  z3  Xj) 


Thus,  the  perturbation  in  CT  is  also  linearly  proportional  to  s. 

An  alternative  method  of  determining  SP  location  As  is  to  solve  for  the  minimum  of  the 
integrand  of  equation  (46)  on  the  real  A-axis.  Equivalently,  the  minimum  of  equation  (47), 
which  is  the  logarithm  of  the  integrand,  will  suffice.  For  the  latter  function,  there  follows 

A  (Xu)  =  A(AS ,  u)  +  A '(/I,  ,u)  (A-As)  +  \A*(As,u)  (A  -  As  )2 
for  A  near  As.  But,  use  of  the  relations 

,w)  =  -w  +  x'(As )  =  0,  A”(AS ,  u)  =  z”(As ), 
leads  to  the  result 


A  (A.  u)  =  A(A,  ,u)  +  \  x”(As )  (A  -  As  )2 .  ( 

Therefore,  when  a  relative  error  tolerance  is  imposed  on  the  search  for  the  minimum,  namely. 


A(/U)-A(;s,t/) 
A  (As,u) 


<R„. 


the  requirement  yields 
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(65) 


I  i-K 


^2^jA(^oiy/2 

,  zW  J 


Now,  let  solution  Ac  (the  axis  crossing  actually  used  for  the  Bromwich  contour  C)  be 


A  -  A  +  e, 


where 


M < 

(2Rt  A(a„«)|'j 

l  zW  J 

(66) 


Then,  the  analysis  in  equations  (53)  through  (56)  is  unaltered,  except  that  now 


R2 


£  z"M 
2  lU) 


|  A  (A’") 


x\kY J 


1/2 


(67) 


is  the  relative  error  of  the  estimated  SPA  pc(u).  A  comparison  with  equation  (57)  reveals  that 
relative  error  R\  depends  linearly  on  tolerance  Ra,  whereas  relative  error  Ri  varies  as  the  square 
root  of  tolerance  Rt-  This  degradation  in  the  quality  of  Rj  is  due  to  the  fact  that  searching  for  a 
minimum  of  a  function  is  fundamentally  less  accurate  than  searching  for  a  zero  of  its  derivative. 
The  minimum  search  is  at  the  bottom  of  a  parabolic  bowl,  where  much  larger  changes  in  the 
abscissa  are  required  than  at  the  zero  crossing  of  a  nearly  linear  function. 

The  bounds  in  equations  (57)  and  (67)  depend  on  high-order  derivatives  of  %(A)  =  log //(A) 

evaluated  at  the  SP  location  As.  All  of  these  quantities  depend  upon  the  value  of  the  field  point  u 
as  well  as  the  particular  PDF  p(u )  under  investigation.  No  general  statements  about  accuracy  are 
apparent;  accordingly,  a  couple  of  examples  will  be  investigated  quantitatively. 


EXAMPLE  1:  EXPONENTIAL  RANDOM  VARIABLE 

p(u)  =  exp(-w)  for  u  >  0, 

p(A)  =  — for  Re(A)  <  1, 

1  —  A 

A  (A,  u)  =  -A  u  -  log(l  -  A), 

/(2)  =  -log(l-i),  z'(A)  =  T—r,  =  for »  >  1 . 

1  —  A  (1  —  A) 


(68) 
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All  these  relations  require  Re(/t)  <  1 .  For  u  >  0,  the  SP  equation  yields 
=  r^-r  =  w»  X=\--, 


^U)  =  7T7^  =  w2’  xmW=2u\  rw  =  6u4, 

V~~As) 

Z’U,)  2 

raf  »' 

Then,  the  relative  error  in  equation  (57)  is  upper-bounded  by 

r,<rA-=i 


That  is,  the  relative  error  in  estimated  SPApc(w)  is  less  than  that  in  determining  SP  location  Xs 
itself,  by  means  of  the  SP  equation  z'(Xs)  =  u. 

By  use  of  equations  (68)  and  (69),  there  follows 


(»-!)! 


X„  =  X  W  =  = 

Zl_2lL  =  24^_22^  = 

X2  6«4  u2 

2  - 3^1  =  2^- - 3—  =  0 

*3  *2  2m3  U2 

Thus,  the  correction  term  (60)  is  unaffected;  that  is, 

1 _  A  =  I  _  A  =  3  _  5  _  _\_ 

8  xW  24  x"U)3  8  1/(1  -X)A  24  1/(1 -A)6"  4  6  _  12 

for  all  X.  The  correction  term  is  independent  of  X  or  Xs  or  Xc.  The  SPA  is 


ps{u)  =  exp(-w)  -=  =  1 .0844  exp(-w)  for  u  >  0. 

■Jin 

For  the  alternative  approach,  where  A (X.u)  is  minimized,  the  relevant  relations 
A (X,.u)  =  -Xs  u  -  log(l  -  Xs )  =  1  -  u  +  log(?/) 
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(74) 


(75) 


;T0U2_(2»3)2 
1U)3  O'2)3  ’ 

giving 

R2<(2Rb\l-u  +  log(w)|  )1/2 . 

One  bad  effect  in  this  approach  is  the  dependence  on  7?^ 2  instead  of  7?*;  this  is  due  to  the 

quadratic  behavior  of  A(A,u)  near  SP  location  As.  The  second  bad  effect  is  the  particular 
dependence  on  u,  namely, 


(2Rb  |log(w)|)1/2  asw-^0- 


(2  Rbuf 


as  u  -»  +oo 


Although  the  singularity  near  u  =  0  is  not  likely  to  cause  problems  because  of  its  weak  behavior, 
the  singularity  for  large  u  could  be  problematical,  especially  if  u  is  the  order  of  1/7?*  »  1.  Even 
if  u  =  0.01/7?*,  the  relative  error  is  (2/100)1/2  =  0.14,  which  is  a  noticeable  relative  error. 


EXAMPLE  2:  GAMMA  RANDOM  VARIABLE 


«'*&-«)  for„>0j 
r(& +i) 

=  forRe(A)<l, 

zW  =  ~(k  +  l)log(l-A), 

A(7l,  u)  =  -X  u  -  (k  + 1)  log(l  -  A), 

j'(A)  =  — ,  j(w>(A)  =  (il  +  l)  for  «  >  1 . 

1-A  (1-/1)" 


For  w  >  0, 


^  + 1  15  ^  +  i  5  £+1 

z(2s)  =  u,  - — r"w’  l~As= — >  K  =1 - 1 

1  -  w  u 

T  (*  +  l)2  ’  J  (  4  ~  (*  +  l)3 

■  ,  independent  of  A. 

/U)2  * 
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Then,  the  relative  error  in  equation  (57)  is  upper-bounded  according  to 

independent  of  k. 

At  the  same  time,  by  the  use  of  equations  (78)  and  (79), 


xn=x(n)(A)  = 


(»-1)!kw 

ot+ir1 


giving  rise  to 


— -2— =  0,  2— -3— =  0. 


X  4  Xi 


Xi  Xi 


Thus,  the  correction  term  is  unaffected;  that  is, 

1  z"W  5  zW .  1  1 

8  zvy  24  zW  12k  +  l 

for  all  A.  The  correction  term  is  independent  of  A  or  As  or  Ac. 

The  SPA  is 

uk  exp(-w)  exp(k  +  1)  T(k  +  1)  ,  .  „ 

Ps («)  =  — 7±=  tTT)  =  P(u)  8(k)  for  u  >  0. 

r(*  +  l)  y/2x  (k  +  \)k*]  2 

In  particular, 

g(0)  =  1.0844,  g(l)  =  1.0422,  g{2)  =  1.0281,  g(10)  =  1.0076,  g(oc)  =  1, 

.  1  23  3821 

g(*)  -  1  +  — —  -  + - -  as£— >+oo. 

\2k  288A2  51840  it3 

For  the  alternative  approach,  where  A (A,w)  is  minimized. 


A(/?,,w)  =  -As  u-(k  +  1)  log(l  -  AJ-  -u  +  k  +  \-(k  +  \)  log 


rUj  _  4 

/U)3  *  +  i 
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which  yields 


R2  < 


2Rh 


1  — 


u 


k  +  1 


+  log 


u 


£  +  1 


,1/2 


(88) 


In  particular, 


R , 


2Rh 


u 


,1/2 


k  +  l 


as  u  — +  +oo. 


(89) 


That  is,  larger  k  values  of  the  Gamma  PDF  ameliorate  the  accuracy  problems  present  for  large 
values  of  u.  The  Rb2  dependence  is  still  present,  however. 


It  has  been  observed  numerically  that  the  SPA  is  very  accurate  for  evaluating  the  joint  PDF 
of  interest  here;  in  fact,  the  correction  term  is  virtually  1  in  most  cases.  The  reason  for  this  is  that 
the  integral  in  equation  (3)  pertains  to  the  last  RV,  which  is  the  sum  r  of  the  remainder.  For  a 
large  number  N  of  RVs,  and  for  M  not  too  large,  the  Central  Limit  Theorem  says  that  r  (the  sum 
of  N  +  1  -  M  RVs)  tends  to  a  Gaussian  RV.  But  the  SPA  itself  is  exact  for  Gaussian  RVs. 
Therefore,  since  r  is  nearly  Gaussian,  it  is  expected  that  the  SPA  should  be  quite  accurate;  this  is 
indeed  the  case. 
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PROGRAMS  FOR  THE  SADDLEPOINT  APPROXIMATIONS 


Six  MATLAB  programs  for  the  SPAs  in  equations  (8)  and  (9)  are  listed  below.  They  are 
keyed  to  the  four  examples  listed  in  the  Saddlepoint  Approximation  section  and  to  tables  1 
through  6. 


Example  1 : 

Example  2: 

Example  3: 

Example  4  (chi-squared,  four  DOF): 
Example  4  (chi-squared,  six  DOF): 
Example  4  (chi-squared,  eight  DOF): 


spa_M_log_l  (table  1) 
spa_M_log_2  (table  2) 
spa  M  log_3  (table  3) 
spa_M_log_4_cs4  (table  4) 
spa_M_log_4_cs6  (table  5) 
spa_M_log_4_cs8  (table  6). 


The  log  reference  in  the  file  names  indicates  that  the  logarithm  of  the  SPA  is  actually  calculated, 
rather  than  the  SPA  itself.  This  avoids  underflow  and  overflow  when  extreme  tail  probabilities 
are  of  interest.  The  abbreviation  csp  stands  for  chi-squared  of  p  DOF,  p  integer. 

The  logarithm  of  the  standard  SPA  is  given  by  variable  spaO  log,  while  the  corresponding 
first-order  correction  term  is  given  by  spal  log.  The  error  factor  was  defined  just  below 
equation  (57).  Fourth-order  derivatives  of  the  various  l  functions  are  required  to  evaluate  these 
expressions.  These  quantities  are  thoroughly  considered  in  appendixes  A,  B,  and  C. 
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Table  1.  Program  for  Example  1 


function  spa_M_log_l 
global  zM  zMl  bd 

N=1024 ; 
n=[l:N] f ; 

b-exp (  - .  3* (n-.7*N) .  A2)+.l; 


%  TR  11,469,  general  M 
%  Exponential  RVs,  CS  2  DOF 

%  Number  of  RVs 

%  N  means 


m=  [10; 11; 12; 13] ;  z=[1.3; .8; .6; .4; 100] ; 


M=length (z) ; 
zd=dif f  (z)  ; 

if (any (zd (1 :M-2) >0) )  pdf=0,  keyboard,  end 
zM=z (M) ; 
zMl-z  (M-l)  ; 

if (zM> (N+l-M) *zMl)  pdf-0,  keyboard,  end 

nd=setdiff (n,m)  ;  %  n  m(l) , . . . ,m(M-l) 

bd=b (nd) ; 


options=optimset ( 1 tolx 1 , le-15)  ; 

[Ls,  fv, ex, ou] =fzero (@Laml, 0, options) ; 
saddlepoint=Ls 

ap=z (1 :M-1)  . /b (m)  ;  %  Leading 

pm_log=-sum (ap) -sum (log (b (m) ) ) ;  %  PDFs 

as=zMl . /bd~zMl*Ls ; 

Lam=-zM*Ls+sum  (log  ( zMl .  /bd)  )  -f  sum  (Lcs2z  (as)  )  ; 
Lam2=zMl A2*sum (Lcs2z2 (as) ) ; 

Lam3=-zMlA3*sum (Lcs2z3 (as)  )  ; 

Lam4=zMl A4*sum (Lcs2z4 (as)  )  ; 
spa0_log=pm_log+Lam- . 5*log (2*pi*Lam2 ) 
ct=.125*Lam4/Lam2A2-5/24*Lam3A2/Lam2A3; 
spal_log=spa0_log+ log ( 1+ct ) 
error_factor=. 5*zM*abs (Lam3) /Lam2A2 
keyboard 

function  w  =  Laml  (L)  %  Lambda1  (L) 

global  zM  zMl  bd 

w=-zM-zMl*sum (Lcs2zl ( zMl . /bd-zMl*L) ) ; 
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Table  2.  Program  for  Example  2 


function  spa_M_log_2  %  TR  11,469,  general  M 

global  zM  zMl  bd  %  Gaussian-squared  RVs 

N=128;  %  Number  of  RVs 

n= [1 :N] ' ; 

b=exp (- . 3* (n- . 7 *N) . A2 ) + . 1 ;  %  n  scale  factors,  page  A-3 

m= [10; 11; 12 ; 13] ;  z= [1 . 3; . 8; . 6; . 4; 10]  ; 

M=length (z) ; 
zd=diff (z) ; 

if (any (zd (1 :M-2 ) >0) )  pdf=0,  keyboard,  end 
zM=z (M) ; 
zMl=z (M-l) ; 

if (zM> (N+l-M) *zMl)  pdf=0,  keyboard,  end 

nd=setdiff (n,m) ;  %  n  ~=  m(l) , . . . ,m(M-l) 

bd=b  (nd) ; 

options=optimset ( 'tolx* , le-10)  ; 

[Ls, fv, ex, ou] =f zero (@Laml , 0, options) ; 
saddlepoint=Ls 

ap=.5*z ( 1 ; M-l ) . /b (m) ;  %  Leading  PDFs 

pm_log=-sum(ap)-.5*sum(log(2*pi*z(l:M-l) .  *b  (m) ) ) ; 

as=. 5*zMl . /bd-zMl*Ls; 

Lam=-zM*Ls+. 5*sum(log ( . 5*zMl . /bd) ) +sum(Lerfz (as) ) ; 

Lam2=zMl A2*sum(Lerfz2 (as) ) ; 

Lam3=-zMl A3*sum (Lerf z3 (as) ) ; 

Lam4=zMl A4*sum(Lerfz4 (as) ) ; 
spaO_log=pm_log+Lam- . 5*log (2*pi*Lam2 ) 
ct=. 125*Lam4/Lam2A2-5/24*Lam3A2/Lam2A3; 
spal_log=spaO_log+log (1+ct) 
error_factor=. 5*zM*abs (Lam3) /Lam2A2 
keyboard 


function  w  =  Laml (L)  %  Lambda' (L) 

global  zM  zMl  bd 

w=-zM-zMl * sum (Lerf zl ( . 5*zMl . /bd-zMl*L) ) ; 


Table  3.  Program  for  Example  3 


function  spa_M_log_3 
global  zM  zMl  ad  bd 

N=128 ; 
n= [ 1 : N ] 1 ; 
a=l+n/N; 

b=exp (  - .  3* (n-.7*N) . A2)+.l; 


%  TR  11,469,  general  M 
%  Gaussian  RVs 

%  Number  of  RVs 

%  N  shift  factors,  page  A- 3 
%  N  scale  factors,  page  A-3 


m=  [117; 126; 89] ;  z= [2 . 5; 2 . 1;2 . 05; 185] ; 


M=length (z) ; 
zd=dif f (z)  ; 

if ( any ( zd ( 1 : M-2 ) >0 )  )  pdf : =0 
zM=z  (M)  ; 
zMl=z  (M-l)  ; 

if (zM> (N+l-M) *zMl )  pdf=0. 


,  keyboard,  end 
keyboard,  end 


nd=setdiff (n,m)  ;  %  n  ~=  m(l)  ,  .  .  .  ,m(M-l) 

ad=a (nd) ;  bd=b (nd) ; 


options=optimset ( 1 tolx1 , le-10) ; 

[Ls, fv, ex, ou] =fzero (@Laml, 0, options) ; 
saddlepoint=Ls 

ap= (z (1 :M-1) -a (m) ) . /b (m) ;  %  Leading  PDFs 

pm_log=- . 5*sum(ap. A2) -sum (log (sqrt (2*pi) *b (m)  ) ) ; 

As=(zMl-ad) . /bd-bd*Ls;  As2=As.A2; 
phis= . 5*erf core (-As* sqrt ( . 5)  ,  1)  ; 
rs=sqrt ( . 5/pi) *exp (- . 5*As2 ) . /phis; 
bd2=bd. A2;  bd2r=bd2.*rs; 

Lam=-zM*Ls+sum(Ls*ad+. 5*LsA2*bd2+log (phis) ) ; 

Lam2=sum(bd2-bd2r .  *  (As+rs)  )  ; 

Lam3=sum(bd. *bd2r . * (l-As2-rs . * (3*As+2*rs)  )  )  ; 

Lam4=sum(bd2 . *bd2r . * (As . * (3-As2 ) +rs . * ( 4-7*As2-rs . * ( 12*As+6*rs) ) ) ) ; 

spaO_log=pm_log+Lam- . 5*log (2*pi*Lam2 ) 

ct= . 125*Lam4/Lam2  A2~5/2  4  *Lam3 A2/Lam2 A3 ; 

spal_log=spaO_log+log ( 1+ct ) 

error_factor= . 5*zM*abs (Lam3) /Lam2 A2 

keyboard 


function  w  =  Laml (L)  %  Lambda*  (L) 

global  zM  zMl  ad  bd 

A= (zMl-ad) . /bd-bd*L; 

phiA= . 5*erf core (-A*sqrt ( . 5) , 1 ) ; 

r=sqrt ( . 5 /pi ) *exp (- . 5*A. A2 ) . /phiA; 

w=-zM+sum (ad+bd . A2*L-bd . *r ) ; 
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Table  4.  Program  for  Example  4,  Four  DOF 


function  spa_M_log_4_cs4 
global  zM  zMl  bd 

N=1024 ; 
n= [1 :N] 1 ; 

b=exp(-.3* (n-.7*N) .A2)+.l; 


%  TR  11,469,  general  M 
%  Chi-squared  RVs,  4  DOF 

%  Number  of  RVs 

%  N  means 


m=  [10; 11; 12; 13] ;  z= [1 . 3; . 8 ; . 6; . 4 ; 170]  ; 

M=length  (z) ; 
zd=dif f (z) ; 

if (any(zd(l:M-2) >0) )  pdf=0,  keyboard,  end 
zM=z (M)  ; 
zMl=z (M-l); 

if (zM> (N+l-M) *zMl )  pdf=0,  keyboard,  end 


nd=setdiff (n,m) ;  %  n  ~=  m(l) , . . . ,m(M-l) 

bd=b (nd) ; 


options=optimset ( ' tolx ’ , le-15) ; 

[Ls, fv,ex,ou] =f zero (@Laml, 0, opt ions) ; 
saddlepoint=Ls 

ap=z ( 1 : M-l ) . /b (m) ;  %  Leading  PDFs 

pm_log=-sum (ap) -sum (log (b (m) ) ) +sum(log (ap) ) ; 

as=zMl . /bd-zMl*Ls; 

Lam=-zM*Ls+2*sum  (log (zMl . /bd) ) +sum(Lcs4z (as) ) ; 
Lam2=zMlA2*suit)(Lcs4z2  (as)  )  ; 

Lam3=-zMl A3*suin (Lcs4z3  (as)  )  ; 

Lam4=zMlA4*sum(Lcs4z4 (as) ) ; 
spaO_log=pm_log+Lam- . 5*log (2*pi*Lam2) 
ct= . 1 2 5 * Lam4 / Lam2 A 2 - 5 /2 4  * Lam3 A 2 /Lam2 A 3 ; 
spal_log=spaO_log+log (1+ct) 
error_factor=.5*zM*abs (Lam3) /Lam2A2 
keyboard 

function  w  =  Laml  (L)  %  Lambda'  (L) 

global  zM  zMl  bd 

w=- zM- zMl *sum ( Lcs  4  z 1 (zMl . /bd-zMl*L) ) ; 


Table  5.  Program  for  Example  4,  Six  DOF 


function  spa_M_log_4_cs6 
global  zM  zMl  bd 

N=1024 ; 
n= [1 :N] 1 ; 

b=exp (  - .  3* (n-.7*N) . A2)+.l; 


%  TR  11, 469,  general  M 
%  Chi-squared  RVs,  6  DOF 

%  Number  of  RVs 

%  N  means 


m=  [10; 11; 12; 13] ;  z= [1 . 3; . 8; . 6; .4; 230] ; 


M=length (z) ; 
zd=diff (z); 

if (any (zd(l:M-2) >0) )  pdf=0,  keyboard,  end 
zM=z  (M)  ; 
zMl=z  (M-l)  ; 

if (zM> (N+l-M) *zMl)  pdf=0,  keyboard,  end 

nd=setdiff (n,m) ;  %  n  m (1) , . . . ,m(M-l) 

bd=b (nd) ; 

options=optimset ( * tolx 1 , le-15) ; 

[Ls, fv, ex, ou] =f zero (@Laml, 0, options) ; 
saddlepoint=Ls 

ap=z (1 :M-1)  . /b (m)  ;  %  Leading  PDFs 

pm_log=-sum(ap) -sum(log (b (m)  ) ) +2*sum(log (ap) ) - (M-l) *log (2)  ; 

as=zMl . /bd-zMl*Ls; 

Lam=-zM*Ls+3*sum(log  (zMl .  /bd)  )  +sum(Lcs6z  (as)  )  ; 
Lam2=zMlA2*sum(Lcs6z2 (as) ) ; 

Lam3=-zMlA3*sum(Lcs6z3 (as) ) ; 

Lam4=zMlA4*sum(Lcs6z4  (as)  )  ; 
spaO_log=pm_log+Lam-  .  5*log  (2*pi*Lam 2 ) 
ct= .  125*Lam4/Lam2 A2-5/24*Lam3A2/Lam2A3; 
spal_log=spaO_log+log (1+ct) 
error_factor= . 5*zM*abs (Lam3) /Lam2  A2 
keyboard 

function  w  =  Laml (L)  %  Lambda*  (L) 

global  zM  zMl  bd 

w=-zM-zMl*sum (Lcs6zl ( zMl . /bd-zMl*L) ) ; 
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Table  6.  Program  for  Example  4,  Eight  DOF 


function  spa_M_log_4_cs8  %  TR  11,469,  general  M 

global  zM  zMl  bd  %  Chi-squared  RVs,  8  DOF 

N=1024;  %  Number  of  RVs 

n= [1 :N] ' ; 

b=exp (-.3* (n-. 7*N) . A2) +. 1;  %  N  means 

m= [10; 11; 12 ; 13] ;  z= [1 . 3; . 8; . 6; . 4 ;270] ; 

M=length (z) ; 
zd=diff (z) ; 

if ( any ( zd ( 1 : M-2 ) >0 ) )  pdf=0,  keyboard,  end 
zM=z (M)  ; 
zMl=z (M-l) ; 

if (zM> (N+l-M) *zMl)  pdf=0,  keyboard,  end 

nd=setdiff (n,m) ;  %  n  ~=  m(l) , . . . ,m(M-l) 

bd=b (nd) ; 

options=optimset ( 'tolx' ,  le-15)  ; 

[Ls,  fv, ex, ou]=f zero (@Laml,0, options) ; 
saddlepoint=Ls 

ap=z ( 1 : M-l ) . /b (m) ;  %  Leading  PDFs 

pm_log=-sum(ap) -sum (log (b (m)  ) ) +3*sum(log (ap) ) - (M-l) *log (6) ; 

as=zMl . /bd-zMl*Ls; 

Lam=-zM*Ls+4*sum(log (zMl . /bd) ) +sum(Lcs8z (as) ) ; 
Lam2=zMlA2*sum(Lcs8z2 (as) ) ; 

Lam3=-zMlA3*sum(Lcs8z3 (as) ) ; 

Lam4=zMlA4*sum(Lcs8z4 (as) ) ; 
spaO_log=pm_log+Lam- . 5*log (2*pi*Lam2 ) 
ct=. 125*Lam4/Lam2A2-5/24*Lam3A2/Lam2A3; 
spal_log=spaO_log+log (1+ct) 
error_factor=. 5*zM*abs (Lam3) /Lam2A2 
keyboard 

function  w  =  Laml (L)  %  Lambda' (L) 

global  zM  zMl  bd 

w=-zM-zMl*sum(Lcs8zl (zMl . /bd-zMl*L) ) ; 


SUMMARY 


Equations  have  been  derived  for  the  standard  and  first-order  corrected  SPA  to  the  combined 
probability  and  joint  PDF  of  the  M—  1  largest  RVs  of  a  set  of  N  independent  RVs  and  the  sum  of 
the  remainder.  These  results  have  been  applied  to  six  different  examples  of  RVs,  including 
Exponential,  Gaussian-Squared,  Gaussian,  and  Chi-Squared  RVs  of  various  DOF.  Detailed 
investigations  of  the  necessary  auxiliary  functions  and  their  singularities  have  been  conducted, 
and  highly  accurate  programs  have  been  generated  for  their  efficient  calculation.  Most  of  the 
results  have  13  to  14  decimal  digits  of  accuracy. 

It  has  been  observed  numerically  that  the  SPA  is  very  accurate  for  evaluating  the  joint  PDF 
of  interest  here;  in  fact,  the  first-order  correction  term  is  virtually  1  in  most  cases.  The  reason  for 
this  is  that  the  integral  in  equation  (3)  pertains  to  the  last  RV,  which  is  the  sum  r  of  the 
remainder.  For  a  large  number  N  of  RVs,  and  for  Mnot  too  large,  the  Central  Limit  Theorem 
says  that  r  (the  sum  of  N  +  1  -  M  RVs)  tends  to  a  Gaussian  RV.  But  the  SPA  itself  is  exact  for 
Gaussian  RVs.  Therefore,  since  r  is  nearly  Gaussian,  it  is  expected  that  the  SPA  should  be  quite 
accurate.  This  is  indeed  the  case. 
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APPENDIX  A 

ALTERNATIVE  POWER  SERIES  FOR  log[(l  -  exp(-z))/z] 


The  function  L(z )  of  interest  and  its  first  four  derivatives  were  given  in  equations  (12)  and 
(16),  namely. 


1(2)  =  log 


f  l-exp(-z) 


for  z*  0,  X(0)  =  0, 


(A-l) 


and 


L\z) 

L\z) 

L\z) 

L'\z) 


+  E  =  exp(z), 


z2  (E-l)2  ’ 

2  E(E  + 1) 
z3  +  (E  - 1)3  ’ 

6  E(E2+4E  +  1) 

z4  (E  -l)4 


(A-2) 


where  z  is  complex.  The  values  of  the  functions  in  equation  (A-2),  at  z  —  0,  are  to  be  taken  as  the 
limits  as  z  approaches  zero.  All  of  these  equations  have  numerical  problems  as  z  approaches 
zero,  and  the  problem  becomes  most  severe  for  the  highest-order  derivative.  The  required 
differences  of  large  numbers,  of  approximately  equal  size,  cause  losses  of  significance  in  the 
final  results  near  z  =  0,  if  equations  (A-l)  and  (A-2)  are  used  in  their  current  forms. 

An  alternative  approach  is  to  replace  these  equations  with  their  power  series  expansions 
about  z  =  0.  Questions  then  arise  as  to  how  many  terms  should  be  used  in  each  power  series,  and 
at  what  values  of  |z|  the  switches  should  be  made,  from  the  given  functional  forms  above,  to  the 
appropriate  power  series.  The  answers  to  these  questions  hinge  on  the  desired  accuracy,  both 
absolute  and  relative,  and  the  amount  of  computer  effort  to  be  expended  on  the  power  series 
equivalents. 

The  approach  to  this  problem  will  be  illustrated  with  respect  to  function  L(z)  in  equation 
(A-l).  First,  “exact”  values  of  L(z)  were  obtained  directly  from  equation  (A-l)  (at  any  nonzero  z 
of  interest)  by  means  of  variable  precision  arithmetic  using  50  digits  of  significance.  Then,  a 
(nonzero)  value  of  radius  r  =  \z\  was  specified  and  L(z)  was  evaluated  via  equation  (A-l)  at  128 
equally  spaced  points  on  a  circle  of  radius  r  in  the  complex  z-plane,  using  standard  double 
precision  (64  bits).  The  maximum  absolute  and  relative  errors  of  these  latter  calculations 
(relative  to  the  “exact”  50-digit  results)  were  then  evaluated  and  plotted  as  one  point  on  each  of 
the  curves  labeled  “Z,”  in  figure  A-l .  The  radius  r  was  then  varied  from  0.05  to  1 .2  in  increments 
of  0.05  and  the  complete  “I”  curves,  both  absolute  and  relative  in  figure  A-l,  were  swept  out. 
These  error  curves  basically  decrease  with  radius  r  because  the  distance  from  the  problematical 
point  z  =  0  is  increasing. 


A-l 


Next,  tiie  power  series  expansion  of  I(z)  to  8  terms,  namely  up  to  z8,  was  determined.  For 
radius  r,  this  power  series  in  z  was  then  evaluated  in  standard  double  precision  at  the  same  128 
points  as  above,  and  the  maximum  absolute  and  relative  errors  (relative  to  “exact”)  were 
computed  and  plotted  in  figure  A-l  with  the  label  “z5”.  As  radius  r  was  increased,  the  errors 
basically  increased  because  the  power  series  is  most  accurate  near  the  origin  of  the  complex  z- 
Plane.  This  procedure  was  repeated  for  the  power  series  expansions  using  10, 12, 14,  and  16 
terms,  and  the  corresponding  absolute  and  relative  errors  are  indicated  in  figure  A-l. 

From  the  plots  in  figure  A-l ,  it  is  now  possible  to  determine  where  to  make  the  switch  from 
the  functional  form  in  equation  (A-l)  to  one  of  the  various  power  series  indicated.  For  present 
purposes  with  function  L(z)  in  equation  (A-l),  it  was  decided  to  make  the  switch  at  radius  r  = 
0.56  and  to  use  the  power  series  with  terms  up  to  z12  for  |z|  <  0.56.  Figure  A-l  indicates  that  the 
relative  error  is  then  better  than  10‘15  for  all  complex  z.  A  larger  switch-radius  would  require  a 
larger  power  series,  but  the  improvement  in  relative  error  would  be  minimal. 

Similar  results  for  the  four  derivative  functions  in  equation  (A-2)  are  presented  in  figures 
A-2  through  A-5,  respectively.  The  quantity  L{p\z)  is  denoted  by  Lp  in  the  figures.  It  was 
decided  to  use  13  terms  for  £(z)  below  radius  0.63, 16  terms  for  L\z)  below  radius  0.9, 23 
terms  for  Lm(z)  below  radius  1.42,  and  24  terms  for  L'\z)  below  radius  1.45.  The  errors  are 
getting  progressively  poorer  for  the  higher-order  derivative  functions.  However,  the  worst 
relative  error  for  L"\z)  still  maintains  about  12.7  decimals,  according  to  figure  A-5.  This 
bound  is  applicable  for  all  complex  z. 

The  slow  decays  of  the  I-curves  with  radius  r  discourage  the  use  of  higher-order  series 
expansions.  However,  if  larger  errors  can  be  tolerated,  figures  A-l  through  A-5  indicate  the 
fr^de-offs  available  to  a  user.  For  example,  if  20  terms  are  acceptable  for  the  power  series  for 
L"\z) ,  the  relative  error  would  be  increased  to  about  12.4  decimals  and  the  switch  radius  would 
be  changed  to  1.2. 

All  of  these  functions  have  poles  at  z  =  ilm  for  n  -  ±1,  ±  2,  ....  This  non-analytic  behavior 
means  that  the  power  series  expansions  about  the  origin  z  =  0  cannot  possibly  converge  for 
[z| >  because  that  is  the  location  of  the  poles  closest  to  the  origin.  Furthermore,  one  cannot 

get  too  close  to  |z|  =  2 n  because  the  convergence  would  then  be  so  slow  as  to  be  futile.  Also 

observe  that  the  pole  of  L""(z)  at  z  —  ±/2;r  is  of fourth  order;  that  singularity  makes  for  very 
slow  convergence  if  |z|  is  allowed  to  get  anywhere  near  2 n. 

Programs  for  all  the  functions  in  equations  (A-l)  and  (A-2)  are  presented  in  tables  A-l 
through  A-5.  These  MATLAB  programs  are  vectorized  for  any  size  matrix  r  and  can  handle  all 
complex  values  of  z. 
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Table  A- 1.  Program  for  L(z)  in  Equation  (A-l) 


function  w  =  Lcs2z(z)  %  log ( [1-exp (-z) ] /z) 
s=size(z);  z=z ( : ) ;  %  for  all  complex  z 

w=zeros (size (z) ) ;  %  and  all  size(z). 

%  Principal  value  log  is  desired.  The  branch  lines  are  at: 

%  z  =  lambertw (n, exp (1/p) /p) -1/p  for  real  p  >  0;  n=+-l,  +-2,  . 
%  Branch  points  of  Log(...)  are  at  z  =  i  2  pi  n,  n=+-l,  +-2,  . 

kl=find(real (z)<=-37) ; 
if (~isempty (kl) ) 
y=z (kl) ; 

t  y-log (-y) ;  %  Not  necessarily  the  principal  value  for 

b=imag(t); 

b=mod (b, 2*pi) ; 

k=find (b>pi) ; 

b(k)=b(k)-2*pi; 

w(kl) =real (t) +b*i;  %  Principal  value  for  w. 

end 

k2=find (abs (z) <=. 56) ;  %  Radius  0.56 

if (~isempty (k2) ) 

y=z(k2);  x=y.*y; 

al=-l/2;  a2=l/24 ;  a4=-l/2880;  a6=l/181440;  a8=-l/9676800; 
al0=2 . 087675698786810e-9;  al2=-4 . 4 034 917 822 3 9578 e-11; 
w (k2) =al*y+x . * (a2+x. * (a4+x. * (a6+x . * (a8+x . * (al0+x.*al2) ) ) ) ) 

end 

k3=setdiff ( [l:length(z) ] ' ,union (kl, k2) ) ; 
if (~isempty (k3) ) 
y=z (k3) ; 

w(k3)=log( (1-exp (— y) ) . /y) ; 

end 

w=reshape (w, s) ; 

_ 


A 


Table  A-2.  Program  for  Lr(z)  in  Equation  (A-2) 


function  w  =  Lcs2zl(z)  %  -1/z+l/ (exp (z) -1) 
s=size(z);  z=z ( : ) ;  %  for  all  complex  z 

w=zeros (size (z) ) ;  %  and  all  size(z). 

%  Function  has  poles  at  i  2  pi  n  for  n=+-l,  +-2,... 

kl=find (real (z) >=41)  ; 
if (-isempty (kl)  ) 

w(kl)=-l./z(kl)  ; 

end 

k2=find (abs (z) <=. 63) ;  %  Radius  0.63 

if (-isempty (k2) ) 

y=z ( k2 ) ;  x=y.*y; 

a0=-l/2;  al=l/12;  a3=-l/720;  a5=l/30240;  a7=-l/1209600; 
a9=l/ 47900160 ;  all=-5 . 284190138 6874 93e-10; 
al3=l . 3382536530 684 68 e-11; 

w(k2) =a0+y. * (al+x. * (a3+x. * (a5+x. * (a7+x. * . . . 

(a9+x. * (all+x. *al3)  ))))); 

end 

k3=setdif f ( [ 1 : length ( z ) ] 1 , union ( kl , k2 ) ) ; 
if (~isempty (k3) ) 
y=z (k3) ; 

w (k3 ) =-l . / y+1 . / (exp(y)-l) ; 

end 

w=reshape  (w,  s) ; 


Table  A-3.  Program for  L"(z)  in  Equation  (A-2) 


function  w  =  Lcs2z2(z)  %  l/zA2-E/ (E-l ) A2;  E=exp(z) 
s=si-ze(z);  z=z  ( : )  ;  %  for  all  complex  z 

w=zeros (size  (z) ) ;  %  and  all  size(z). 

%  Function  has  poles  at  i  2  pi  n  for  n=+-l,  +-2, . . . 

kl=find (real (z) >=45) ; 
if (~isempty (kl) ) 

w(kl)-l./z(kl) .  A2; 

end 

k2=find (abs (z) <=. 9) ;  %  Radius  0.9 

if {~isempty (k2) ) 

y=z(k2);  x=y . *y; 

a0=l/12;  a2=-l/240;  a4=l/6048;  a6=-l/172800;  a8=l/5322240; 
al0=-691/118879488000;  al2=l/5748019200; 
a!4=- 3617/71137 48561 92 000;  al6=1.459630549567234e-13; 
w(k2 ) =a0+x. * (a2+x . * (a4+x . * (a6+x . * (a8+x. * . . . 

(alO+x.* (al2+x . *  (al4+x.*al6) )))))); 

end 

k3=setdif f ( [1; length (z) ] ' , union (kl, k2) ) ; 
if (~isempty (k3) ) 

y=z(k3);  E=exp(y); 
w(k3)=l./y.A2-E./(E-l) .A2; 

end 


w=reshape (w, s) ; 


Table  A-4.  Program  for  L'"(z)  in  Equation  (A-2) 


function  w  =  Lcs2z3(z)  % 
s=size(z);  z=z ( : ) ;  % 
w=zeros (size  (z) ) ;  % 
%  Function  has  poles  at  i 


-2/zA3+E (E+l) / (E~l) A3;  E=exp(z) 
for  all  complex  z 
and  all  size (z) . 

2  pi  n  for  n=+-l,  +-2,... 


kl=find (real (z) >=49) ; 
if (~isempty (kl) ) 

w(kl)=-2./z (kl) . A3 ; 

end 


k2=f ind (abs (z) <=1 . 42 ) ;  %  Radius  1.42 

if (~isempty (k2) ) 

y=z(k2);  x=y.*y; 

al=-l/120;  a3=l/1512 ;  a5=-l/28800;  a7=l/665280; 

a9=-691/11887948800;  all=l/47 9001600;  al3=-7 . 118328 622277424e-ll; 
al5=2 .335408879307 57 4 e- 12;  al7=-7.438050949068572e-14; 
al9=2. 313781187 911296e-15;  a21=~7 . 060959131021136e-17 ; 
a23=2 . 12 0824223777 681e-18; 

w(k2) =y. * (al+x. * (a3+x. * (a5+x. * (a7+x. * (a9+x. * (all . . . 

+x.* (al3+x. * (al5+x . * (al7+x. * (al9+x . * (a21+x . *a23) )))))))))); 

end 


k3=setdif f ( [1 : length (z ) ] f , union (kl, k2 ) )  ; 
if (-isempty (k3) ) 

y=z ( k3 ) ;  E=exp (y) ; 
w(k3)=-2./y.A3+E.*(E+l) . /(E-l) .  A3; 

end 

w=reshape  (w,  s)  ; 
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Table  A-S.  Program  for  L""(z)  in  Equation  (A-2) 


function  w  =  Lcs2z4 (z)  %  6/zA4-E(EA2+4E+l) / (E-l) A4;  E=exp(z) 

s=size (z) ;  z=z ( : ) ;  %  for  all  complex  z 

w=zeros (size (z) ) ;  %  and  all  size(z). 

%  Function  has  poles  at  i  2  pi  n  for  n  =  +-1, 

kl=f ind ( real ( z ) >=52 ) ; 
if (-isempty (kl) ) 

w(kl)=6./z(kl) . A4 ; 

end 

k2=f ind (abs (z) <=1 . 45) ;  %  Radius  1.45 

if (-isempty (k2) ) 

y=z(k2);  x=y.*y; 

a0=-l/120;  a2=l/504;  a4=-l/5760;  a6=l/95040; 
a8=-691/1320883200;  al0=l/43545600;  al2=-9. 253827208 960651e-10; 
al4=3 . 503113318961361e-ll;  al6=-l .264468661341657e-12; 
al8=4 . 3961842570314 63e-14;  a20=-1.482801417514439e-15; 
a22=4 . 8778 95714 688 665e-17;  a24=-1.571342308445089e-18; 
w (k2 ) =a0+x . * (a2+x. *  (a4+x . * (a6+x.* (a8+x.* (alO+x.* (al2. . . 

+x. * (al4+x. * (al6+x. * (al8+x. * (a20+x. * (a22+x. *a24) )))))))))); 

end 

k3=setdif f ( [1 : length (z) ] ' , union (kl, k2) ) ; 
if (~isempty (k3) ) 

y=z(k3);  E=exp(y); 

w(k3)=6./y.A4-E.*(E.A2+4*E+l) ./(E-l) .A4; 

end 

w=reshape (w, s) ; 
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APPENDIX  B 

ALTERNATIVE  POWER  SERIES  FOR  logferf(Vz  )/vz] 


The  function  L(z)  of  interest  and  its  first  four  derivatives  were  given  in  equations  (19)  and 
(23),  namely. 


L{z)  =  log 


for  z  *  0,  1(0)  =  log(2/ 4n), 


(B-l) 


where  z  is  complex.  The  values  of  the  functions  in  equation  (B-2)  at  z  =  0  are  to  be  taken  as  the 
limits  as  z  approaches  zero.  All  of  these  equations  have  numerical  problems  as  z  approaches 
zero,  and  the  problem  becomes  most  severe  for  the  highest-order  derivative.  The  required 
differences  of  large  numbers,  of  approximately  equal  size,  causes  loss  of  significance  in  the  final 
results  near  z  =  0,  if  equations  (B-l)  and  (B-2)  are  used  in  their  current  forms. 

The  method  of  circumventing  this  numerical  difficulty  is  identical  to  that  presented  in 
appendix  A,  namely,  replace  these  functions  by  their  power  series  representations  for  z  near  0. 
The  results  for  the  absolute  and  relative  errors  for  the  five  functions  above  are  presented  in 
figures  B-l  through  B-5,  respectively.  For  the  /.-function  error  in  figure  B-l,  it  was  decided  to 
use  equation  (B-l),  as  is,  for  all  z.  On  the  other  hand,  it  was  decided  to  use  6  terms  in  the  power 
series  forZ/(z)  below  radius  0.06,  16  terms  for  L”(z)  below  radius  0.71,  18  terms  for  L”(z) 
below  radius  0.87,  and  22  terms  for  L"\z)  below  radius  1.08.  Thus,  for  example,  from  figure 
B-5,  the  relative  error  for  L"”(z)  is  12.1  decimals  for  all  complex  z. 


There  are  bumps  in  some  of  the  relative  error  curves.  These  are  due  to  the  exact  function 
having  a  zero  for  a  particular  value  of  z.  For  example,  Z-(z)  has  a  zero  at  z  =  0.3812.  Lm(z)  has  a 
zero  at  z  =  - 1 .4585,  and  L"\z )  has  a  zero  at  z  =  0.9087.  When  the  relative  error  is  computed  at 
a  point  near  these  zeros,  these  near-zero  calculations  of  the  denominator  manifest  themselves  as 
bumps  in  the  relative  error  curve,  and  should  be  ignored. 
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MATLAB  programs  for  all  five  functions  in  equations  (B-l)  and  (B-2)  are  available  from 
the  author  upon  request. 
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Figure  B-5.  Errors  for  L""(z)  in  Equation  ( B-2 ) 
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APPENDIX  C 

ALTERNATIVE  POWER  SERIES  FOR  CHI-SQUARED  RANDOM  VARIABLES 


For  k  =  1,  that  is,  for  a  chi-squared  RV  with  2k  +  2  =  4  DOF,  the  function  L(z)  and  its  first 
four  derivatives  were  given  in  equation  (38),  namely. 


m = log 


1  -  exp(-z)  (1  +  z) 


forz*0,  Z(0)  =  -log(2). 


Z'(z)  =  — -  +  — ^ ,  E  =  exp(z), 

z  E-l-z 

r(2)=4+Mz£M 

z2  (E-l-z)2 

L”(z)  =  -±+  —  (~2  +  2)  +  -£(4-z  +  z2)-2 
^  *  z3  (E-l-z)3 

t ffff/  \  12  £3(3-z)  +  4£2(-3  +  2z-z2)  +  £(15-7z  +  z2-z3)-6 

L  (Z)=—  H - ; - . 

z4  (E-l-z)4 


(C-l) 


where  z  is  complex.  The  values  of  the  functions  in  equation  (C-l),  at  z  =  0,  are  to  be  taken  as  the 
limits  as  z  approaches  zero.  All  of  these  equations  have  numerical  problems  as  z  approaches 
zero,  and  the  problem  becomes  most  severe  for  the  highest-order  derivative.  The  required 
differences  of  large  numbers,  of  approximately  equal  size,  causes  losses  of  significance  in  the 
final  results  near  z  =  0,  if  the  results  in  equation  (C-l)  are  used  in  their  current  forms. 

The  method  of  circumventing  this  numerical  difficulty  is  identical  to  that  presented  in 
appendix  A,  namely,  replace  these  functions  by  their  power  series  representations  for  z  near  0. 
The  results  for  the  absolute  and  relative  errors,  for  the  five  functions  above,  are  presented  in 
figures  C-l  through  C-5,  respectively.  It  was  decided  to  use  13  terms  in  the  power  series  for  L(z) 
below  radius  0.7, 14  terms  for  L'(z)  below  radius  0.9,  15  terms  for  L"(z)  below  radius  1.03, 18 
terms  for  Lm(z )  below  radius  1.27,  and  23  terms  for  L""{z)  below  radius  1.84.  Thus,  for 
example,  from  figure  C-5,  the  maximum  relative  error  for  L""(z )  is  12  decimals  for  all  complex 


The  relative  error  curve  for  L(z)  has  a  bump  in  the  neighborhood  of  radius  1 ;  this  is  due  to  a 
zero  of  L(z)  at  z  =  -1 .  Also,  Lm(z )  has  a  zero  at  z  =  2. 1 63 1 ,  while  L""(z)  has  zeros  at 
z  =  - 1 . 1 756  and  5.4429;  however,  the  only  one  of  these  latter  zeros  that  affects  the  relative  error 
curve  is  the  bump  in  the  neighborhood  of  radius  1 .2  in  the  plot  of  L"\z ) . 
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For  *  ~  2’  *at  is,  for  a  chi-squared  RV  with  2k  +  2  =  6  DOF,  the  function  L(z )  and  its  first 
four  derivatives  were  given  in  equation  (39),  namely 


L{z)  =  log 


l-exp(-z)  (1  +  2  +  z2/ 2) 


for  z  *  0,  Z(0)  =  -  log(6), 


3  z2 


Z'(z)-  ,  E  s  exp(z),  D  =  E-(l  +  z  +  z2/2), 

^  ,  E(2  —  z)  —  2  —  z 

iw-7+z — W - ’ 


6  ,  £2(4-8z  +  2z2)  +  £(-8  +  8z  +  8z2-2z3  +  z4)  +  4-6z2-2z3 

z3  4Z)3 

=  ~  +  ^-(4£3(-6  +  6z-z2)  +  8Z2(9-8z2+4z3-z4) 

+  £(-72  -  72z  +  68z2  +16z3  -14z4  +2z5  -z6)  +  6(4  +  8z-4z3  -z4)). 


(C-2) 


The  numerical  difficulties  at  z  =  0  are  handled  as  discussed  above;  the  results  for  the  absolute 
and  relative  errors  for  these  five  functions  are  presented  in  figures  C-6  through  C-10, 
respectively.  It  was  decided  to  use  13  terms  in  the  power  series  for  L(z)  below  radius  0.91, 13 
terms  for  L\z)  below  radius  1.06,  15  terms  for  L\z)  below  radius  1.33, 23  terms  for  Z*(z) 
below  radius  2.15,  and  25  terms  for  Ln\z)  below  radius  2.38.  Thus,  for  example,  from  figure 
C-10,  the  maximum  relative  error  for  V\z)  is  1 1.5  decimals  for  all  complex  z. 


The  relative  error  curve  for  L(z)  has  a  bump  in  the  neighborhood  of  radius  1 ;  this  is  due  to  a 
zero  of  L(z)  at  z  =  -0.9289.  Also,  Lm(z)  has  a  zero  at  z  =  3.9977,  while  L’\z)  has  a  zero  at 

z  =  7.6874;  however,  the  plots  did  not  have  to  be  carried  out  that  far  to  make  the  choices  for  the 
switch  radius. 
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For  k  =  3,  that  is,  for  a  chi-squared  RV  with  2k  +  2  =  8  DOF,  the  function  L(z)  and  its  first 
four  derivatives  were  given  in  equation  (40),  namely 

.  .  fl-exp(-z)(l  +  z  +  z2/2  +  z3/6)^l  f  r/m  i 

L{z)  =  log  - — — — — - - —  for  z*0,  L( 0)  =  -log(24), 

l  2  J 

4  7^ 

L'(z)  = - + - ,  £  =  exp(z),  D  =  £-(l  +  z  +  z2/2  +  z3/6), 

z  6D 

L\z)  =  4"f  (£  (6  “  2z^  "  (6  +  4z  +  z2 )), 

z  12  D 

Lm{z)  =  -\  +  -^-T  (< 6E 2  (6  -  6z  +  z2)  +  £(-72  +  24z2  +  12z3  -  3z4  +  z5 )  (C-3) 

z3  36  D3  v 

+  36  +  36z  +  6z2  -  12z3  -  6z4  -  z5 ), 

L”\z)  =  — - — —  (36£3  (6-18z  +  9z2  -z3)  +  24E2  (-27  +  54z  +  27 z2  -  15z4  +  6z5  - z6) 

216D4  v 

+  E (648  -  648z  -  1620z2  -1188z3  +126z4  +162z5  +3z6  -21z7  +3z8  -z9) 

- 3(72-216z2  -336z3  -186z4  -24z5  +20z6  +8z7  +z8))+^. 

z 

The  numerical  difficulties  at  z  =  0  are  handled  as  discussed  above;  the  results  for  the  absolute 
and  relative  errors  for  these  five  functions  are  presented  in  figures  C-l  1  through  C-15, 
respectively.  It  was  decided  to  use  15  terms  in  the  power  series  for  L(z )  below  radius  1.45,  16 
terms  for  L\z)  below  radius  1.66,  17  terms.for  L\z)  below  radius  1.88,  22  terms  for  Lm(z) 
below  radius  2.58,  and  23  terms  for  L"”(z)  below  radius  2.90.  Thus,  for  example,  from  figure 
C-15,  the  maximum  relative  error  for  V”(z)  is  10.9  decimals  for  all  complex  z. 

MATLAB  programs  for  all  the  functions  in  equations  (C-l)  through  (C-3)  are  available  from 
the  author  upon  request. 
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